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Chapter 1 


Introduction 


In forging, the material is squeezed between two or more dies in such a way 
that the metal or alloy is plastically deformed to the desired shape and size. 
It is also used to improve the metallurgical and mechanical properties of the 
material. This is the oldest of the metal working processes known to mankind 
since the copper age. 

Generally the forging is performed by heating the material and then 
applying a compressive force to obtain the final shape. This is called the 
hot forging. The cold forging operation is carried out at room temperature. 
The forging operation is carried out in tw'o different ways : drawing out and 
upsetting. 

There are four types of forging methods w'hich are generally used. 

Smith forging- This is the traditional forging operation done openlv or 
in open dies by a black smith, by manual hammering or by power hammers. 
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Dtop forging - This is the operation done in closed impression dies by 
means of the drop hammers. 

Press forging- Similar to drop forging, the press forging is also done 
in closed impression dies with the exception that the force is a continuous 
squeezing type applied by the hydraulic presses. 

Machine forging- Unlike the drop or press forging where the material is 
drawm out, in machine forging, the material is only upset to get the desired 
shape. 

A majority of the finished products made by forging are geometrically 
complex and the metal flow invoked is of a three-dimensional nature. Thus, 
any technique of analysis will become more useful for problems of industrial 
importance, if it is capable of solving three dimensional metal flow problems. 

Residual stresses are self-equilibrating internal stresses existing in a body 
after the removal of external forces and are caused by non-uniform deformation. 
It is generally believed that poor shape or bad surface finish of the work- 
piece is caused by residual stresses. Further, a tensile residual stress field 
near the surface cein cause a micro-crack to propagate rapidly. Experimental 
methods for determination of residual stresses have been developed and used 
[1]. How'ever, w'hen the process is to be optimized, experimental methods do 
not provide an economical solution. Therefore, it is necessary' to determine 
them by analytical or numerical methods. To calculate residual stresses in 
forging, a large deformation elasto-plastic finite element analysis, based on the 
updated Lagrangian formulation, needs to be carried out. 
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1.1 Literature Survey 


The problem of upsetting of metals has theoretical as well as technological sig- 
nificance in that solutions would help to bring about a better understanding of 
the behavior of metals during plastic deformation at high pressure. Such solu- 
tions w'ould give specific answers to problems actually encountered in practice. 
It is for these reasons that a number of investigators have dealt with the forging 
problem. Early investigators considered the material as an ideal plastic solid. 
This leads to so-called slip line solution. Such solutions have been proposed 
by Prandtl [2], Hill et al. [3], and Green [4]. The slip line theory has been 
well developed to analyse a non-homogeneous plane strain deformation of a 
rigid-perfectly plastic isotropic solid. This theory can be used to get very good 
first-approximations to loads required for performing operations and provide 
indications of the manner in which material deforms. However, the elastic 
and work hardening effects can not be incorporated in this method. There- 
fore a detailed solution for the deformation that takes place during continued 
loading beyond the yield point can not be obtained using this technique. Fur- 
ther, residual stresses also cannot be studied using this technique. Shabaik [5] 
investigated the bulging in upsetting by using the slipline theory. 

Hoffman and Sach [6] proposed the slab method. A complete analysis of 
the slab method has been presented by Altan [7] for the axis 3 Tnmetric closed 
die forging. In this method, a slab of infinitesimal thickness is selected and 
the equation representing the force equilibrium for the slab is derived, assum- 
ing that the deformation is homogeneous within the slab and that the major 
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directions are the principal stress directions. The resulting differential equa- 
tion is then solved with appropriate boundary conditions. While this method 
can give very good predictions of the load variation with deformation, it is 
inherently incapable of predicting shape changes such as barreling in open die 
forging. Further, since the method is primarily bcised on the assumption of 
homogeneous deformation, it possibly cannot predict the residual stresses as 
they are caused only by non-uniform deformation. 

The upper bound theorem was formulated by Prager and Hodge [8]. If 
surfaces of velocity discontinuities are included [9], then it states that among 
all kinematically admissible velocity fields r*, the actual one minimizes the 
following expression 

r = I^S*j€*jdV + [ r\Av'\ds- [ TiV^ds (1.1) 

Here is the strain rate field derived from u’, S-j is the deviatoric 
stress field derived from and | Av* \ is the discontinuity in the tangen- 
tial component of the velocity across the surface Sr- Further r denotes the 
shear stress on the surface Sr and Ti is the prescribed traction on the part 
St of the boundary. The first term expresses the pow'er spent in causing the 
internal deformation over the volume V of deforming body. The second term 
represents the power loss over the surfaces of velocity discontinuities including 
the boundary between tool and material. The last term represents the power 
supplied by specified tractions. What the upper bound theorem states is that 
the actual externally supplied power is never higher than that computed by 
the first two terms of the above equation. 
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Kudo [10] applied the upper bound theorem to the problems of axisym- 
metric cold forging and extrusion. Kudo [11] also applied the upper bcund 
theorem to the analyses of the plane strain forging and extrusion. Avitzur 
[12], Kobayashi [13], McDonald et al. [14] and many others have suggested up- 
per bound velocity fields to predict the forming load in upsetting with bulging. 
Yang and Kim [15], Kim et al. [16] and Marques and Martins [17] proposed up- 
per bound method for analysis of three dimensional upset forging. The major 
difficulty in using the upper bound method is how to choose the kinematically 
admissible velocity field since the accuracy of the solution depends on how- 
close the assumed velocity field is to the actual one. Another disadvantage of 
this method is that the solution provided by the method does not satisfy the 
equilibrium equations. So, even if elastic effects are included in integral (1.1), 
the method possibly can not be used for determination of residual stresses. 

Accurate determination of forging parameters under realistic conditions 
became possible when the finite element method was introduced. The major 
advantage of the finite element method is that the method can be applied 
to a wide class of boundary value problems w-ithout restrictions on workpiece 
geometrj- used in the finite element method. Using FEM, it is possible to 
predict the platen forces, the interface pressure and the stress, strain and 
residual stress distributions within the billet at \-arious deformation levels. 

Early applications of the FEM to forging problems were based on the 
incremental method proposed by Lee and Kobayashi [18]. The method uses 
the elastic-plastic stress strain matrix based on Prandtl-Reuss equations. The 
additivity of incremental elastic and plastic strains is assumed. Even though 
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the stress-strain matrix and the geometry are updated after every increment, 
only the linearized incremental equations are used. The method was applied to 
solid cylinder upsetting [19], ring compression [20] and for predicting defects 
in upsetting [21]. Maccarini et al. [22] used the method for studying the 
influence of die geometry on cold extrusion forging operation. Where as Lee 
and Kobayashi [18] used the velocity as the primary unknown, Hartley et 
al. [23-24] proposed an incremental method with the displacement as the 
primary unknown. In their method also, the linearized incremental equations 
are used and the elastic-plastic matrix and the geometry are updated after 
every increment. They studied the upset forging process by incorporating the 
friction using the beta-stiffness method [24].' Three dimensional elasto plastic- 
finite element analysis has been done by Pillinger et al. [25] using linearized 
incremental equations. 

Linearized incremental equations give only an approximate solution. If 
the incremental size is not sufficiently small, the error between the exact and 
the approximate solutions grows rapidly with the applied load, as the error in 
the solution of the current increment gets propagated into the next increment 
when the stress-strain matrix and the geometry are updated. To awid this 
phenomenon, the non-linear incremental equations should be solved by using 
a scheme like Newdon-Raphson technique. Such a formulation was first pro- 
posed by Bathe et al. [26]. This formulation ( with or without elastic effects ) 
has been applied to the forging problem by many researchers. It was applied 
by Carter and Lee [27] to axisymmetric upsetting. The three dimensional 
rigid- viscoplastic finite element analysis was done by Park and Kobayashi [28]. 
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Recently Choi et al. [29] developed three dimensional rigid-plastic finite ele- 
ment formulation with special attention to treatment of the contact between 
the dies and workpiece. At high reduction, a part of the free surface of the 
workpiece comes in contact with die. This phenomenon has been included 
in the formulation of Choi et al. [29]. Terziyski et al. [30] used Arbitrary' 
Lagrangian-Eulerian (ALE) FEM for three dimensional forging. Surdon and 
Chenot [31] developed rigid viscoplastic finite element formulation for hot forg- 
ing process. It considered thermal effects but didn’t use objective stress and 
strain measures. It neglected elastic effects as the rigid plastic/viscoplastic ap- 
proach converges in fewer iterations as compared to elastic plastic/viscoplastic 
approach. Fu and Luo [32] used rigid-viscoplastic finite element formulation 
for predicting defects in isothermal forging. 

The thrust of most of the above studies has been to find an accurate 
estimation of the forging load, to study the deformation field or to predict 
fracture. As far as the determination of residual stresses is concerned, there 
are hardly any attempts [33]. A major difficulty in finding the residual stresses 
in forging is that the increment size has to be very small in order to achieve a 
reasonable level of accuracy in the stress values. This leads to a prohibif -ely 
large amount of computational time. The Jaumann stress measure, which is 
the objective stress measure used presently in the literature, does not give 
accurate stress values if the incremental shear deformation is large. Thus, 
there is a need to look for a new objective stress measure which wilt allow the 
use of a large increment size without compromising on the accuracy of stresses. 
Such a measure has been recently proposed in reference [34]. 
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1.2 Modelling of Cold Forging Process 

In the present work, residual stresses in three dimensional cold upset forging 
are studied. The updated Lagrangian formulation which is convenient for ’can- 
dling material and geometric nonlinearities is used. New incremental objective 
stress measure and logarithmic strain measure are used instead of Jaumann 
stress rate and Green-Lagrange nonlinear strain used in most of the literature. 
Forging process is a displacement controlled problem. Therefore, to accelerate 
the convergence of iterative scheme, arc length method is used in conjunc- 
tion with modified Newton Raphson iterative technique. During unloading, 
the process can be modelled as a force controlled problem, hence only modi- 
fied Newton Rapson iterative scheme is used. The material is assumed to be 
elastic-plastic strain hardening yielding according to von Mises criterion. The 
effect of temperature and strain rate (viscoplastic effects) on the yield strength 
of the material are ignored in this work. The inclusion of these effects renders 
the analysis quite complex. Due to small acceleration, inertial forces are not 
included. The body forces are also neglected. 

The dies are assumed to be rigid and rough. So no relative displacement 
is permitted at the die-workpiece interface and the process is considered un- 
der condition of complete sticking at the die-workpiece interface. The sliding 
friction at the die-workpiece interface, which applies to the case of smooth 
lubricated die, is not considered due to lack of time. 
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1.3 Objective 


The objective of the present work is to develop a 3-D large deformation eleisto- 
plastic finite element code for determination of residual stresses in cold forging 
process. The updated Lagrangian formulation is to be used along with new 
incremental objective stress measure and logarithmic strain measure. The 
arc length method is to be used in conjunction with the Newton-Raphson 
technique to solve the non-linear incremental equations of the displacement 
controlled problem. 

First, the code is checked for mesh convergence and then validated by 
comparing the results with experimental results of reference [28]. Next, a 
detailed parametric study of residual stresses is carried out to show the ef- 
fects of three process variables namely the percentage reduction in height, the 
geometric ratios ( height to width ratio, width to thickness ratio ) and the 
material properties. For a typical case, the normal stress distribution at the 
die-w'orkpiece interface and the deformation patterns after loading and unload- 
ing are also presented. 

1.4 Plan of The Thesis 

The thesis is organized as follows. 

In the second chapter, the mathematical modelling of forging process, 
the 3-D finite element formulation and the numerical scheme are presented. 
The boundary conditions are discussed at the end of chapter. - 
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In chapter 3, the mesh convergence of the code, the \'alidation and the 
results for a typical case are presented. It also includes a discussion on the 
parametric study. 

Conclusions and suggestions for further work are presented in chapter 4 . 
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Chapter 2 


Mathematical Modelling and 
Finite Element Formulation of 
3-D Forging Process 


In this chapter mathematical model and Finite Element formulation of a cold 
forging process is Developed. The process is modelled as a 3-D problem. The 
description of motion, stress measures, strain measures and constitutive ’ iws 
used in the formulation of the governing equation for static large deformation 
elasto-plastic problems are presented in sections 2. 1-2.6 . The finite element 
formulation and numerical schemes, required to implement the finite element 
formulation, are presented in sections 2. 7-2. 9. The boundary conditions for 
3-D forging are discussed in section 2.10. 

2.1 Updated Lagrangian Formulation 

In the study of the deformation of a body subjected to external loading, oden 
the original undeformed and unstressed state of the body is used for the for- 
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mulation of its equation of motion. This is known as Lagrangian formulation. 
This formulation is convenient for small deformation which is the case in a 
majority of engineering problems. In such cases, the deformed configuration 
does not deviate much from the original one and hence the deformation can 
be described by an infinitesimal strain tensor, known as the engineering strain 
tensor, for which the strain-displacement relations are linear. On the other 
hand, for large deformation problems, one has to use a finite strain measure 
which is expressed by a nonlinear strain-displacement relation. Furthermore, 
the equations of motion when expressed in the reference configuration depend 
on the deformation. Hence, for all large deformation problems, the Lagrangian 
formulation proves to be cumbersome with the governing equations being dif- 
ficult to solve. In such cases, one solves the problem using an incremental 
method known as Updated Lagrangian formulation. In this formulation, it is 
assumed that the state of the stress and deformation of the body is knowm till 
the current configuration, say at time t. The main objective is to determine 
the incremental deformation and stresses during an infinitesimal time step, 
At, i.e. from time t to t -F At. Here, the current configuration is used as the 
reference configuration for obtaining the incremental veilues. Unlike in the La- 
grangian formulation, an incremental strain tensor is used. This methodology 
is particularly useful for the elasto-plastic materials because the stress-strain 
relationship in such materials is usually expressed in an incremental fashion. 
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2.2 Kinematics of Finite Deformation 

The Cartesian co-ordinate ij of a particle at time r with reference to the 
co-ordinate Xj at time t is given by 

x, = h(Xi,T) ( 2 . 1 ) 


w'here t < t < t + At. 

The relative deformation gradient at time t + At is defined by = 

fij and if is not singular, the polar decomposition theorem [35] allow's 

a decomposition of the form: 

( 2 . 2 ) 

w'here is an orthogonal tensor representing the material rotation -and 

is a positive definite symmetric tensor called the right stretch tensor. 
The right stretch tensor can be diagonalized by the following transformation 
to obtain the principal stretches : 

= tQik tQjl t^kl 

where is orthogonal and 

(no sum) (2.4) 

The tensor and the transformation matrix will be used h» the 

development of a stress measure while the incremental logarithmic strain def- 
inition will be obtained from the 
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2.3 Stress Measures 


It is essential that an objective stress measure be used in the incremental theory 
of constitutive modeling to account for the rigid body translation and rotation 
that may accompany deformation. The Cauchy stress tensor which has great 
physical significance is not objective and hence cannot be used directly in 
a constitutive equation. There are numerous objective stress measures and 
algorithmic formulations (See [36], [37] and the references therein), each with 
particular advantages and disadvantages. Some of them are discussed below. 


The second Piola-Kirchoff stress tensor S can be related to the Cauchy 
stress tensor a using the concept of equivalent work between two configurations 
[38]: 


t+At c ^ 

tOij — 


t-¥At^i,rn 


t^At 


^mn t+At^ji^ 


(2.5) 


where x represents position vector and p denotes density. The Piola-Kirchoff 
stress tensor is energy conjugate with the Green-Lagrange strain tensor. This 
energ)^ conjugate pair is used in the equation for predicting displacemeni,s in 
a predictor-corrector solution procedure (see section 2.7). 


Jaumann discussed the formulation of the rate of change of the stress 
tensor. The J-rate a is related to the Cauchy rate & by- 

‘ Uii At = ^aij At - At) - At) (2.6) 

where Ml'y At = | {tAuij —t Au^,,) represents the components of the incre- 
mental spin tensor. Here, tAu is the incremental displacement vector. 

As pointed out by Metzger and Dubey [39], it is important' that the stress 
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measure used be compatible with the constitutive equation in addition to being 
objective. An incremental objective stress measure (as against objective ^ess 
rate measures) to be used in the generalized Hooke’s law and in the elasto- 
plastic constitutive equation is developed below. 



The fixed frame is denoted by subscript F 
The material frame is denoted by subscript M 

Figure 2.1: Fixed and material reference frames 

Two Cartesian reference frames (See Figure 2.1) are used : 

1. The fixed frame. 

2. A material frame which rotates and translates along with a material 
particle. 

The material frame is defined so as to coincide with the principal axes of the 
right stretch tensor at time t so that the initially orthogonal axes do not get 
skewed at time t + At. 
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The first step in the derivation is the transformation of the components 
of the Cauchy stress tensor at time t from the fixed frame to the material 
frame: 

(2.7) 

The increment in the components of the Cauchy stress tensor with respect to 
the material axes, (Aa,y, is added to to obtain the stress components at 
time t + At with respect to the material frame of reference : 

= <^ij + tAOi^ (2.8) 

The final transformation from the material frame at time t + At to the fixed 
frame gives the components of the Cauchy stress tensor at time t + At : 

t+At _ t+At/^ t+At p t+At _M t+At p t-fAt/^ /o q\ 

— t^ki t^kl <^lm fVnj K^-^) 

The equations (2.7), (2.8) and (2.9) can be combined together to give a relation 
between the components of the Cauchy stress tensors at times t and t + At 
and the increment in the components with respect to the material frame : 

t+A«_ _ t+Atrt '+AtD ft+Atri t„ t+At^ , i+A£p £+At^ 

^ij — t^ki t^kl ( ^mn t^on "i t^^lo J t^pj 

( 2 . 10 ) 

The proposed stress and the logarithmic strain measures satisfy the re- 
quirements of objectivity and lead to a physically consistent application of the 
usual constitutive equation. 
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2.4 Strain Measures 


The incremental small strain is given by 

. 1 , . 

tAcjj — (jAuj-j + tAzij,) (2.11) 

The incremental Green-Lagrange strain is a non-linear function of the displace- 
ment: 

— (jAujj f -j- (2.12) 

These two strain measures occur in the virtual work expression at time t + At 
and its transformation to time t respectively. The components of the Green- 
Lagrange strain tensor are invariant under rigid body rotation of the material 
unlike the small strain components. 

The following are the disadvantages of using one of the above measures 
in a constitutive law : 

1. The solution obtained is dependent upon the size of the increment in the 
updated Lagrangian formulation unless the increment size is sufficiently 
small. 

2. The components of the strain tensors do not tend to infinite values when 
the principal stretches tend to zero. Therefore a constitutive law which 
ensures that the appropriate Cauchy stress components tend to negative 
infinity (as is physically realistic), even though the strain components 
remain finite, should be used. This difficulty can be avoided by using 
a strain measure whose components become minus infinity wffien the 
principal stretches become zero. 
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The logarithmic strain measure introduced by Dienes [36] is free of the 
above disadvantages. The principal incremental logarithmic strain compo- 
nents, which will be used in this work, are defined by 

= (2.13) 

where the Aj are the principal stretches defined by equation (2.4). The loga- 
rithmic strain has the following additional advantage in elasto-plastic analysis. 
A loading test involving elasto-plastic deformation followed by elastic unload- 
ing reveals that the slope of the elastic unloading line is the same as that of 
the initial elastic loading line only when the true stress and the logarithmic 
strain measures are used in a constitutive law [40]. 

2.5 Elastic Constitutive Equation 

A constitutive law must satisfy the principle of material frame indifference [35]. 
The elastic constitutive equation used is the generalized Hooke’s law relating 
the increment in stress components with respect to the material frame and the 
elastic part ([Ae®-^]) of the principal incremental logarithmic strain components 

,Aa*'=|“'cJ„rf(,A4f) (2.14) 

The tensor for the isotropic case is given by 

X 5,, 5ki+2fi5ik5ji (2.15) 

where A and fi are Lame’s constants. 
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In case anisotropic behavior is modeled, the tensor C® has to be evalu- 
ated with reference to the directions of the principal stretches. 

The stress and strain measures used in a constitutive equation need not 
necessarily be energy conjugate with each other. However if it is so, the pre- 
dicted response in a predictor-corrector scheme will be closer to the actual 
response. 


2.6 Elasto- plastic Constitutive Equation 

As stresses developed in a material exceed the yield stress, the linear relation- 
ship between the stress and strain no longer remains valid. We now develop a 
relationship between the stress and strain based on the von Mises yield crite- 
rion and isotropic hardening ^ . 

For an isotropically hardening material, the plastic potential is given by 

[41] ; 

-f’(t'y.p) = - PjW {2-16) 

Note that 

F = 0 (2.17) 

represents the yield criterion. The plastic potential F depends on the Cauchy 
stress tensor <7^ through it’s second invariant Cgg called as equivalent stress 
and defined by 

(2.18) 

tMs section, the superscript/ subscript denoting time has been omitted for the sake 
of convenience. 
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where is the deviatoric part of a^j. Further, F depends on the variable 
yield stress of the material, ay, through a hardening variable p. For the case of 
strain hardening, p is identified as the equivalent pleistic strain and hence 
defined as : 

= (2.19) 


and 




( 2 . 20 ) 


Here, de^j is the plastic part of the incremental linear strain tensor dcy and 
the integration in equation (2.19) is to be carried along the particle path. The 
dependence of on p (or e^^) is normally approximated by a power-law tj'pe 
of relationship 


a, - (a,)o = K{e%r (2-21) 

Here, ( 0^)0 is the yield stress at zero plastic strain, K is called the hardening 
coefficient and n is called as the hardening exponent. 


The plastic part of incremental linear strain tensor (defj) is obtained firom 
the plastic potential using the following relation : 


dF 

= dX— (2.22) 

^ aoy 

where dX is a scalar. This equation is called as the flow rule. Differentiation 
of equation 2.16 with respect to cfy gives 


dOij 2£Te,^‘^ 

Using this, one can determine dX as : 


(2.23) 


dX = del. 


(2.24) 



Further, the hardening relationship and the yield condition can be used to 


express dX as 


where 


ax- — - — 


H 


day 

del 


KnielX-^ 


(2.25) 


(2.26) 


eg 


is the slope of the hardening curve. Substitution of equations (2.23) and (2.25) 
in equation (2.22) leads to the following constitutive equation 




3 / 

2 0^°“ 


(2,27) 


This constitutive relationship between the deviatoric stress tensor and 
the plastic part of incremental linear strain tensor is not really convenient for 
the Updated Lagrangian formulation for which the incremental stress-strain 
relationship is needed. This can be obtained from equation (2.27) as follows : 


I p 3 a^j da^q 

ddij = - - -^ doki 

^ 2 HOeq OUkl 

Note that, from equations (2.16) and (2.23), we get 

daeg _ dF ^ . 

dau doki 2aeg^’‘‘ 


(2.28) 


(2.29) 


Substitution of equation (2.29) in equation (2.28) leads to the following incre- 
mental plastic stress-strain relationship : 




(2.30) 


The incremental elastic stress strain relationship (equations 2.14-2.15) 
can now' be written as 


de^j = — {-i/dakkSij + (1 + v>)daij] 


(2.31) 
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where defj is the elastic part of dcy, E is the Young’s modulus and u is the 
Poisson’s ratio. Adding the two relationships, w^e get 


diij = de^j + de^j 


r , (l + ^)r r 


eq 


do. 


kl 


(2.32) 

(2.33) 


This is the incremental elasto-plastic stress-strain relationship needed in the 
Updated Lagrangian formulation. How'ever, it is the following inverse relation- 
ship which is more useful : 


doij — C^jj^idc/ci 


(2.34) 


where 


riEP 

'-'ijkl 


2p 


^ik^jl + 


V 


2u 


dijSki — 


9fiaiio‘ 


ij^kl 


(2.35) 


2(3/z + /f)a2,, 

and fi is one of the Lame’s constant ( also called as shear modulus ) appearing 
in eq. (2.15). 


Note that the stress increment appearing in eq. (2.34) must be an objec- 
tive stress increment in the sense that dOij must reduce to a zero tensor in the 
event of the increment being a pure rigid rotation. The incremental objective 
stress measure to be used in the present work has been described in section 
2.3 . 


The relationship (2.34) has been derived assuming the increment size 
to be small and using the incremental linear strain tensor as the strain mea- 
sure. When a large-size increment is to be used along with the incremental 
logarithmic strain measure (defined in section 2.4), the derivation is similar. 
The relationship (2.34), w'hen the stress and strain measures are replaced re- 
spectively by the incremental objective stress measure (of section 2.3) and the 
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incremental logarithmic strain measure (of section 2.4), takes the following 
form with respect to the material frame : 


.A< = I 






(2.36) 


where 


'CfX = 2/x 






l-2u 


^ij^kl 


9 ^ 




2(3/1 + VI,, 


(2.37) 


Here, H has been replaced by the expression (2.26) and the left superscript t 
has been added to make it explicit that these quantities are to be evaluated at 
that time. 


2.7 Incremental Updated Lagrangian Formu- 
lation 

The objective of the updated Lagrangian formulation is to establish static 
equilibrium in the configuration at time t + At when all static variables at 
time t are known. 

The principle of virtual w'ork requires that 

/+Atv (2.38) 

where 

Cartesian component of the Cauchy stress tensor at time t + At 
2 \ ) 

Cartesian component of the virtual displacement -at time t + At 


= 
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— Cartesian co-ordinate of a material point at time t + At 
t+Aty _ Volume at time t -f At 

and 

(2.39)'' 

Jj+A(5^ 

where 

= Cartesian component of external traction per unit area at time t -f At 
t+AtRf — Surface at time t + At with traction specified 

The main difficulty in the application of equation (2.38) is that the con- 
figuration at time t + At is unknown. An elegant way of formulating the 
problem is given by Bathe [42]. The virtual work expression at time t + At 
is transformed to an integral over the volume at time t by using the second 
Piola-Kirchoff stress and Green-Lagrange strain measures, which are energy 
conjugate to each other, and the principle of conservation of mass. It is as- 
sumed that the external load term (2.39) is deformation independent for the 
formulation of the governing equation. The expression (2.38) after the trans- 
formation becomes 

d^V = (2.40) 

The Green-Lagrange strain tensor is defined by 

= \ (2.41) 

The following quantities are incrementally decomposed : 

= JS,,- -h t 

= *<Jij + tASij - (2.42) 
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= 5Cui + tAu.) = SitAui) (2.43) 

= <5(tAe.,) 

= S{tAeij + tArjij) (2.44) 

where 

tA^ij — ((Atijj + lAuj^i) (2.45) 

tAtj-ij = — tAufi^i t Aujtj (2.46) 

Therefore equation (2.40) can be WTitten with incremental decomposition as 

tAS^j SitAtij) d^V + tASij SitArjij) d^V 
+ (5(tAr?y ) d^V + Vy 5{tAtij) d^V = (2.47) 


The above equation is linearized by neglecting the second integral w'hich is a 
higher order term and approximating t A5y as tAtki : 

tA€ki SitAei^) d^V + % SitArjij) d^V 

+ Vij 5(tAey) d'y = (2.48) 

The linearized equation, w^hen solved, will \ield only approximate displace- 
ment, strain and stress fields. The approximate quantities are denoted by a 
right superscript (1). The error due to the approximation involved is calculated 
from equation (2.38) as 

Error = - f 

Jt+Atxr(l) U ' O ^ 

(2.49) 

This error is generally minimized by a predictor-corrector scheme. The next 
section discusses the techniques used for obtaining the solution to the governing 
equation (2.47). 
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2.8 Finite Element Formulation 


A finite element method for solving the linearized equation and an iterative 
numerical technique for minimizing the error is described in this section. 

2.8.1 Matrix Notation 

Matrix notation is used in the development of the finite element equations for 
the convenience of computer implementation. 

The components of the tensors fAe, jAr; and jAc-^ are represented in 
array form as follow's: 

tiAe} = "{{Aej;!, (Acyy, jAe^j, 2i^Cxy, 2tAey;j, 2i^6xz} (2.50) 

~ , jAu^y , f A u^ 2 , fAu^i , i Au^y , • . .}■ (2.51) 

t{A€^} = {tAe^j., tACyy, tAcj^, (2.52) 

The components of all stress tensors are written as arrays w^ith respective 
components placed in the following manner: 

{o’} = {Oli, Oyy, CTzZj O’lyj O’yj , o’^x} (2.53) 

Two matrix forms of the tensor (given by eq. 2.37) are needed: 

1. for the constitutive relationship t{AS} = t{Ac} 

2. for the constitutive relationship t{Ao^} = *(C^'^]t{At^} 
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They are given by 


tf^EP]' _ / rpE'i _ ] *{a} *{c}^[C^ ] \ 

^ ^ ^ ‘H+ ‘{anC^']'{o}/ 

,i^£P. ^ firE, _ jcfmHancf]_\ 

‘ ' V ‘H + ‘{a}^lC^I<{a} J 

In the above equations. 


(2.54) 

(2.55) 


t 


{a} 



(2.56) 


and the array ‘{cr'} represents the deviatoric part of the Cauchy stress at time 
t. For an isotropic material, the matrices [C^ ] and [C^] for the 3-D case are 
given by 



E 

{l + u){l-2u) 


1-1/ 

u 

V 

0 

0 

0 

u 

l-u 

V 

0 

0 

0 

V 

u 

1 - V 

0 

0 

0 

0 

0 

0 

\-2v 

2 

0 

0 

0 

0 

0 

0 

l-2t/ 

2 

0 

0 

0 

0 

0 

0 

1-21- 

2 


(2.57) 


[C^] = 


E 


\ — V 

u 

V 

0 

0 

0 

V 

1-1/ 

V 

0 

0 

0 

u 

V 

1 - V 

0 

0 

0 

0 

0 

0 

l-2y 

0 

0 

0 

0 

0 

0 

l-2t/ 

0 

0 

0 

0 

0 

0 

l-2v 


(2.58) 


Equation (2.48) can be written in the following form owing to the sym- 
metries in jAc, t-^T] and ‘o’ : 


/ i(,{Ae}"')‘[C®'‘'!,{A6}<i‘V+ [ <S(,{A7)}’')‘lr],(A7!})rf‘K 

jty jty 

-I- ^(t{Ae}^) ‘{o} d^V = (2.59) 
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The matrix ^[T] is given by: 


where 




0 

0 ■ 


•1T] = 

0 

‘K 

0 

(2.60) 


0 

0 

‘P], 



CTxx 




•Pl = 



^(^yz 

(2.61) 


^zx 


^ZZ 



2.8.2 Finite Element Equation 


The domain is discretized into a number of elements and the incremental dis- 
placement field is approximated over each element by 


tAu ' 

t{Aw} = •{ tAr 

tAu’ ^ 




(2.62) 


where the element incremental displacement vector t{Au}* is given by 


t{Au}®’' = -[tAu^ tAt;\ tAu!\ tAn^. fAr^, tAtr;^, . . 


(2.63) 


The quantities tAu*, tAu*, tAtu* stand for the unknown incremental displace- 
ments of node i in .iY, Y and Z directions respectively and the matrix *[$] is 
defined by 






(2.64) 


where 


= {‘iVi, 0,0,%. 0,0, 

‘{$2^ = {0,%, 0,0,%.0,0, ...}, (2.65) 

‘{$3}^ = {0,0,%, 0,0,%,0, ...} 


28 



The *Ai, which are functions of (‘ar/y/ z), are called shape functions and are 
determined by the element used. The 8-noded brick element with linear shape 
functions [43] is used in this work. 


The strain field is expressed in terms of the nodal displacements by dif- 
ferentiating (2.62) and using the expressions (2.45) and (2.46) as 


t{Av} = 


where 


and 


w = 


1 

‘{^ 2 } 5 + ml 

+ {^2}, 

ml. 


( 2 . 66 ) 

(2.67) 


( 2 . 68 ) 


= [m^^ 'm,v^ 'm,:^, 'm,^, 'm,v, •--] ( 2 . 69 ) 


Using equations (2.62), (2.66) and (2.67), the contribution to the integral 
(2.59) over a typical element e with volume V‘ is 

S '[Bif ‘[C^’’'] ‘IBil <fV‘),[Av.Y + 

5 (,{A« }'’■)(/_„ '[B«! ^V‘) ,{Au}' + 

S{,{Auf)(l^^ ‘[Btr‘{<T}i‘V'‘ = i{,{Ai.r’')'+'^*{Br (2-70) 

The contribution to the term is expressed in terms of the elemental 

external force vector using a standard procedure [43]. Since the 
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variation in the displacement vector is arbitraiv-, the above equation can be 
written as 


^[K]\{IluY + *{/}' = (2-71) 

The elemental stiffness matrix ^[K]^ is expressed as 

^[Kf = (2.72) 

'\Kl]‘ = (2.73) 

‘[KnlT = f ‘\B^f‘\TY[B„]d!V‘ (2.74) 

The elemental internal force vector is 

‘{/r= / (2.75) 


The elemental stiffness matrices ^[KY along with the elemental force 
vectors ‘{/}^ and are assembled to obtain the global equation ; 

‘[/r]t{Au}+ ‘{/}= (2.76) 

Decomposing t+At{F}, eqn. (2.76) can be written as 

^[K] ,{Au} +‘ {/} = ^F} +, {AF} (2.77) 

Here, ‘{F} is the (global) external force vector at time t and t{AF} is the 
(global) incremental force vector from time t to t+At . In Updated Lagrangian 
formulation, it is assumed that the equilibrium equations are satisfied exactly 
at time t. Thus 

•{/} = '{F} (2.78) 

Then eqn (2.77) reduces to 

^[K]t{Au}= t{AF} ' (2.79) 
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This equation is solved to obtain the incremental displacement vector 
The wavefront equation solver as described by Irons and Ahmad [44], has 
been implemented with buffer and blocking concepts to solve large systems of 
equations. A resolving facility, w'hose need will become apparent later, is also 
present. 

2.8.3 Determination of Stresses 

The evaluation of the stress components (at the Gauss points of the elements) 
is done by the following stepwise procedure : 

1. Calculate the relative deformation gradient 

= 5,, + (2.80) 

It should be noted that the position vector 'x corresponds to the equi- 
librium position. 

2. Decompose and determine and 

using equation (2.3). 

3. Determine the principal incremental logarithmic strain using 

equation (2.13). 

4. Calculate t{AcT^} using equation (2.36). The integration of the consti- 
tutive equation is performed using the Euler forward technique described 
below'. 

5. Use equation (2.10) to calculate the Cauchy stress components at time 

t + At. 
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2.8.4 Integration of the Constitutive Equation 

Different techniques exist for the integration of the constitutive equation [38]. 
A simple and robust technique is the Euler forward integration scheme de- 
scribed below. 

Assume that the principal incremental strain components have been cal- 
culated and the state of the Gauss point at time t (elastic or plastic) is known. 


• If the state at time t is elastic, 


1. Calculate the stress increment assuming elastic behaviour 




2. Calculate using (2.10) 

3. Determine Vg, and using (2.18) 

4. then elastic behaviour holds, or if 
the Gauss point is neutrally loaded. Return. 

5. If a transition from elastic to plastic state has 

occurred. Calculate 


Ratio 




and change the state to plastic. The sub-incrementation method 
is followed and the stress tensor components with ri^pect to the 
material frame are updated for each sub increment by the increment 
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in stress components corresponding to the elasto-plastic strain sub- 
increment. The [C^^] corresponding to last updated state is used: 
t+At{^Ar}W ^ ^ for i = 1, n 

where 

n 

t+Ai{^Ar|( 0 ) ^ + Ratio [C®]t{Ae^} 

t+At^^BP^(o) ^ evaluated at 

Use equation (2.10) to find Return. 

• If the state at time t is plastic, the sub incrementation method described 
in (5) above is applied with Ratio set to zero. 

This describes the Euler forward scheme for the elasto-plastic model described 
in section 2.6 . 

2.8.5 Unloading Scheme 

The phenomenon of local unloading has to be incorporated to reproduce more 
closely the elasto-plastic response of the structure. The unloading criterion 
defined by Chakrabarty [45] is implemented: 

< 0 (2.81) 

where ‘{Ac^} is the incremental strain and ‘{n} represents the unit outwmrd 
normal to the yield surface at the current stress point V in a 9— D space. Neu- 
tral loading and positive loading are represented by replacing the ’’less than” 
symbol in (2.81) by the ’’equal to” and ’’greater than” symbols respectively. 
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When unloading is detected, the yield stress at that Gauss point is 
changed to the equivalent stress at the Gauss point at time t (Vj, = ^o-gg), 
the state tag of the Gauss point is changed from plastic to elastic and the 
elastic constitutive equation is applied to calculate [A<7^^]. 

2.9 Numerical Scheme 

2.9.1 Modified Newton-Raphson Scheme 

The solution to equation (2.79) represents only an approximate equilibrium 
equation at time t + At (the approximation is mostly due to the linearization 
and simplification involved in the steps between eqn (2.47) to (2.48)). A so- 
lution of such an approximate equation may involve a significant amount of 
error and, depending on incremental displacement step, may become unstable. 
Therefore, it is necessary to modify eqn (2.79) to turn it into an iterative prob- 
lem capable of pro\’iding a solution with desirable accuracy. There are various 
iterative techniques for this [38]. The simplest, and for many practical appli- 
cations effective technique is the modified Newton-Raphson algorithm which 
offers fast convergence with less computation. 

This technique can elegantly be represented as follows: 

Solve 

*[A']t{Au}« = (2.82) 

where 

- ‘+^‘{/}b-i) (2.83) 
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(2.84) 


till the convergence criterion 


{F}(0|| 


U+Af 


< toL 


(2.85) 


is satisfied. The vector t{i?} is called as the unbalanced force vector. The right 
superscript on denotes the configuration on which the integration is 

performed to find the external force vector. 


2.9.2 Arc Length Method 


Forging problem is a displacement controlled problem i.e. the deformation is 
controlled by a specified displacement (at the workpiece-die interface) rather 
than by a specified traction. Since there is no knowm external force, the denom- 
inator in the convergence criterion (2.85) doesn’t exist. In that case, one can 
try to achieve the convergence by making the unbalanced force vector t 
small in an absolute sense. But this slows down the rate of convergence to 
a considerable extent. To accelerate the rate of convergence, the arc length 
method is used in conjunction with the modified New^ton-Raphson method. In 
this approach, the unknown nodal reactions at the interface are expressed as a 
linear combination of a set of known nodal forces. Then the problem is solved 
as a force controlled problem using the set of known nodal forces. The un- 
known coefficients in the linear combination are found from the condition that 
the vertical component of the incremental displacement vector at the interface 
has a specified value. Since this solution may not satisfy the nodal equilib- 
rium of internal and external forces, iterations need to be performed which are 
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carried out by the modified Newton-Raphson method. Since both the nodal 
reactions as well as the unbalance force vector change during an iteration, the 
incremental equations contain an additional term compared to eq. (2.82). 

Originally, the arc length method was proposed only for the case of pro- 
portional loading [46, 47]. However, in forging problem, the loading is non- 
proportional. Therefore, the method needs to be modified appropriately to 
take care of non-proportional loading. The modified method is described in 
this section. 

In ith. iteration, the incremental eq. (2.77) is expressed as ^ 

[A^ {Au}^^ = {AFf^ + (2.86) 

where is the unbalanced force vector given by 

= {ry-^ - {/}’-' (2.87) 

and the incremental force vector is expressed as 

m 

{AF}'‘> = y (AA)t [P], (2,88) 

ifc=l 

Here, m is the number of nodes at the interface and {P}k is a known nodal 
force vector corresponding to a unit (vertically downward) point force at node 
k. The vectors {P}^ , k = 1, m are called as the basic load vectors. The 
coefficients (AA)^^ are unknowms. 

' Since eq. (2.86) is a linear equation, one can decompose the solution as 
{Auf^ = f;(AA)f {Au"f^ (2-89) 

A:=:l 

this section, the superscript/ subscript denoting time has been omitted for the sake 
df' convenience. 
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where and ^ are obtained as solutions of the following prob- 

lems : 


[iC] {AuW^} = {P}i, k = l,m 
[K] {Au^^}^‘’ = {RY'-^'> 


(2.90) 

(2.91) 


Thus, represents the iterative displacement vector due to the ba- 
sic load vector at node k and ^ represents the iterative displacement 

vector due to the unbalanced force vector . 


To find the complete solution of eq. (2.86), one needs to determine 
(AA)j\ k = 1, m. They are determined from the following condition. At 
the nodes on the interface, the vertical or z-component of the incremental 
displgicement is known. Thus, for nodes 1 = 1, m : 

m 

^(AA)^'^ = specified value (2.92) 

A:=l 

These are m equations in m unknowns : (AA)[*\ A: = 1, m. By solving these 
equations, the unknown AA^*^ are calculated . Then the iterative displacements 
are determined from eq. (2.89). These displacements are used to find the 
internal force vector {/}^‘^ Finally, the external force vector at the end (of 
ith) iteration is found from 

m 

{F}« = ^ AA« {P}* (2.93) 

k=l 

This force vector is used to find the unbalanced force vector and then to check 
the convergence using eq. (2.85). 

While solving eqs. (2.90-2.91), the upper triangularisation of the coeffi- 
cient matrix [K] needs to be done only once. Then, one can u6e the resohing 
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facility first to perform the Gauss elimination operations on different right side 
vectors and then to find | and ^ by back substitution. 

Equation (2.92) takes a different form for the first iteration than for other 
iterations. These fomrs are as follows : 

First iteration : 

For the first iteration, = {0} implies that ^ = {0}. Thus, 

eq. (2.92) becomes 

m 

^(AA)^’^ = specified value (2.94) 

*=1 

Subsequent iterations : 

Since the value of specified incremental displacement does not change 
during an iteration cycle, in subsequent iterations, the right side of eq. (2.92) 
becomes zero. Thus, 

m 

5](AA)t'^ =0 t > 2 (2.95) 


2.9.3 Numerical Integration Scheme 

Exact evaluation of integrals appearing in element coefficient matrices and 
right side vectors is not always possible because of the complexity of the inte- 
grands. In such cases, it is natural to seek numerical ev’aluation of these integral 
expressions. Numerical integration, involves approximation of the integrand 
by a polynomial of sufficient degree, because the integral of a polynomial can 
be evaluated exactly. 
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The most used numerical integration method is Gauss-Legendre Bumer- 
ical integration method. The expression for the integral over a (3-D) Cubic 
master element is given by 

/n, m 1- 0 d(, dri, = /« /« /+' F(f , C) <i{. dr,, di 

m n I 

- H H Ck)wi ^'k (2.96) 

i=i j=i fc=i 

where m.nj denote the number of quadrature points (Gauss-points) in 77 , C 
directions, and u',. u'j, Wk denote the corresponding weights. 2 Gauss points in 
each direction are used in this work. 

2.9.4 Divergence Procedures 

The modified Newton-Raphson method fails in some cases. Some simple but 
fairly effective procedures for handling divergence are discussed. 

1. Under- relaxation : The iterative displacement is scaled 

<!= 0 <a„ < 1 (2.97) 

where is determined so that the numerical method does not diverge. 

2. Line search : In this technique, the displacement vector is updated ac- 
cording to [38] 

+ a,dAu}^‘^ (2.98) 

and Qj is determined so that the unbalance (2.83) is minimum. 

A full line search is computationally expensive. Therefore a cheaper 
technique is to evaluate unbalance at n different qj values and choose 
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that ai which corresponds to the minimum unbalance. This offers a 
good compromise between computational effort spent in line search and 
effort saved due to convergence in a reduced number of iterations. 

3. Elastic response prediction ; The elastic stiffness is used to predict the 
response and iterations are used to obtain the elasto-plastic response. 
This is beneficial when the structural response is highly stiffening and 
the loading is almost proportional in the increment. 

4. Incremental displacement cutting : In case the procedures described 
above fail, incremental displacement cutting is performed. 



2.10 Boundary Conditions for 3-D Forging of 
Rectangular Block 


Z 

k 



Figure (2.2) shows the domain of the problem. Due to symmetry, only 
1/8 th part of the workpiece is used for analysis. The boundary conditions for 
the domain of this problem are discussed below: 

1 Workpiece-Die Interface (Plane ABCD) 

In this study, it is assumed that the die surface is very rough. So complete 
sticking is assumed at the workpiece-die interface. Then the components of 
incremental displacement in x and y-direction are zero. Since this problem 
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is displacement controlled, the component of incremental displacement in z- 
direction is specified. Thus, boundary conditions on plane ABCD are ; 

Essential boundary conditins : 

dux = 0, duy = 0. du. = specified (2.99) 

2 Free Surfaces (Planes BDFH & CDGH) 

These boundaries are traction free surfaces. Therefore, the boundaiy- 
conditions on planes BDFH & CDGH are : 

Natural boundary conditions : 

dtx = 0. dty = 0. dtz — 0 (2.100) 

3 X-Y Plane of Symmetry (Plane EFGH) 

Because of symmetry, normal component of incremental displacement 
and shear components of incremental traction vector are zero. Thus, the 
boundary conditions on plane EFGH are : 

Natural boundary conditions : 

dtx = 0 . dty = 0 ( 2 . 101 ) 

Essential boundary conditins : 


duz = 0 


( 2 . 102 ) 


4 X-Z Plane of Symmetry (Plane ABEF) 
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Because of symmetry, normal component of incremental displacement and 
shear components of incremental traction vector are zero. Thus, the boundarj' 
conditions on plane ABEF are : 

Natural boundary conditions : 

dt^ = 0 , dt, = 0 ( 2 . 103 ) 

Essential boundary conditins : 

duy = 0 ( 2 . 104 ) 

5 Y-Z Plane of Symmetry (Plzine ACEG) 

Because of symmetry, normal component of incremental displacement 
and shear components of incremental traction vector are zero. Thus, the 
boundary conditions on plane ACEG are ; 

Natural boundary conditions ; 

dty = 0 , dt^ = 0 ( 2 . 105 ) 

Essential boundary conditins : 

du^ = 0 ( 2 . 106 ) 
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Chapter 3 


Results And Discussion 


The finite element model of 3-dimensional forging process developed in the 
previous chapter has been applied to a number of cases involving two materials 
and various sets of input variables to illustrate its applicability. The metals 
considered are Aluminium Al-llOO and Steel AISI 1019 (only in section 3.4.3) 
whose material properties are given as below. The properti^ for the aluminium 
are taken from [28] and that of the steel from [23]. 

Material properties: 

Aluminium Al-llOO 

Young’s modulus E = 69 GPa 

Poisson’s ratio u = 0.3 

Initial yield stress {cry)o = 62.74 MPa 

Hardening coefficient K = 90MPa 
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Hardening exponent n = 0.5377 


Steel AISI 1019 

Young’s modulus E = 210 GPa 
Poisson’s ratio v = 0.3 
Initial yield stress (o-y)o = 478.2 MPa 
Hardening coefficient K = 563.98MPa 
Hardening exponent n = 0.14 

In this chapter, first a convergence study is performed to check the con- 
vergence of the finite element formulation with respect to displacement incre- 
ment and mesh size. Then, the finite element code is validated by comparing 
the results with experimental results. In the last section, the parametric study 
is carried out to illustrate the effects of various process variables on the process 
parameters. 


3.1 Convergence Study 

The convergence study is performed for forging of a block of Al-llOO alu- 
minium. Only 1/8 th of the block is considered for analysis (Fig. 3.1) due 
to symmetries in X, Y and Z planes. First, the convergence of forging load is 
checked. Then the convergence of stress distribution is studied. In each case, 
convergence is checked with respect to refinement in both the mesh size as well 
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Figure 3.1: The domain used for analysis 
as the increment size. 

3.1.1 Convergence of Forging Load 

The block dimensions chosen for this study are : a = b = 9.525 mm and 
h = 4.7625 mm. This size is the same as used in experimental studies of 
reference [28]. The reason for using this size is to facilitate the comparison 
with experimental results of reference [28]. 

Figure 3.2 shows the variation of forging load with % reduction for 4 
different mesh sizes at the increment size of 0.025 mm. When the mesh is 
refined only in z-direction from 5x5x6 elements to 5 x 5 x 12 elements, it 
is observed that there is a large change in forging load (a decrease of 38.5 % 
at 20 % reduction). This happens because the condition of complete sticking 
is assumed at the die-workpiece interface which creates large displacement 
gradients in z direction. When the mesh is refined only in x apd y directions 
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Figure 3.2: Convergence with mesh size 

from 5 X 5 X 12 to 10 X 10 X 12 elements, it is found that there is no change in 
forging load. Thus, the displacement gradients doesn’t seem to be large in x 
or y direction. A further refinement in z direction from 5x5x12 elements to 
5 X 5 X 24 elements brings down the forging load by 14 % at 20 % reduction. 
This study shows that there is a trend towards convergence as the mesh size 
is refined. Analysis with the mesh of 5 x 5 x 24 elements takes a lot of time. 
Therefore, it is decided to use the mesh size of 5 x 5 x 12 elements in the 
further studies. 

Figure 3.3 shows the variation of forging load with % reduction for 3 
different increment sizes for the mesh of size of 5 x 5 x 12 elements. It is 
observed that as the increment size is decreased from 0.05 mm to 0.025 mm, 
the forging load at 20 % reduction decreases by 9.28 %. A further reduction 
of the increment size to 0.0125 mm causes the reduction of forging load (at 
20 % reduction) by 3.7 %. Thus, there is a trend towards convergence as the 
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increment size is decreased. It is decided to use the increment size of 0.025 
mm for the rest of the study. 



Figure 3.3; Convergence with increment size 

3.1.2 Convergence of Stress Distribution 

The block of 9.525 mm x 9.525 mm x 4.7625 mm was used for studying the 
convergence of forging load. The convergence was quite slow. Hence, it was 
decided to reduce the block size for studying the convergence of stress distri- 
bution. So, it was decided to use the size of 5 mm x 5 mm x 5 mm. 

Figures 3.4 to 3.11 show the equivalent stress distribution in x-y and 

y-z planes after 20 % reduction and after unloading for 3 different mesh sizes 

(5x5x6, 5x5x1 2, 10x10x12 elements) and 2 different increment sizes 

(0.025 mm, 0.0125 mm)F It is observed from Figs. 3.4-3.6 that the equivalent 

*The mesh size of 5 x 5 x 6 means 6 elements in z direction and 5 elements each in x and 
f directions. 
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Figure 3.4: Equivalent stress distribution after 20 % reduction in (a) x_y (b) 
y_z planes for 5 x 5 x 6 elements and increment size of 0.025 mm 
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Figure 3.5: Equivalent stress distribution after 20 % reduction in (a) x_y (b) 
V-z planes for 5 x 5 x 12 elements and increment size of 0.025 mm ( I’q 
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Figure 3.6: Equivalent stress distribution after 20 % reduction in (a) x_y (b) 
y_z planes for 10 x 10 x 12 elements and increment size of 0.025 mm ti'n 
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Figure 3.7: Equivalent stress distribution after 20 % reduction in (a) x_y (b) 
y_z planes for 5 x 5 x 12 elements and increment size of 0.0125 mm (fn 

stress values (at 20 % reduction) converge as the mesh size is refined first from 
5x5x6 elements to 5 x 5 x 12 elements and then to 10 x 10 x 12 elements. 
The convergence is also observed when the increment size is decreased from 
0.025 mm to 0.0125 mm (Figs. 3.5 and 3.7). 
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Figure 3.7: Equivalent stress distribution after 20 % reduction in (a) x.y (b) 
y_z planes for 5 x 5 x 12 elements and increment size of 0.0125 mm 

stress values (at 20 7c reduction) converge as the mesh size is refined first from 
5x5x6 elements to 5 x 5 x 12 elenjents and then to 10 x 10 x 12 elements. 
The convergence is also observed when the increment size is decreased from 
0.025 mm to 0.0125 mm (Figs. 3.5 and 3.7). 
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As far as the convergence of residual stresses is concerned, it is observed 
that the pattern of stress distribution in y-z plane changes when the mesh is 
refined in z direction from 5x5x6 elements to 5 x 5 x 12 elements (Figs. 

3.8 and 3.9). The values also decrease drastically. As said in the last section, 
this happens due to large gradients in z direction arising due to the complete 
sticking condition at the interface. A refinement in x and y directions (see Figs. 

3.9 and 3.10) does not cause much change in the pattern. But the maximum 
values do change by as much as 25 %. So the residual stress distribution, even 
though it show's a convergence trend, does not seem to have converged for 
the mesh size of 10 x 10 x 12 elements. It w'as not possible to make further 
refinement. Even the mesh of 10 x 10 x 12 elements takes a lot of time for 
determination of residual stresses at 20 % reduction. So it was decided to use 
the mesh of 5 x 5 x 12 elements with the understanding that the values would 
be off at least by 25 %. When the increment size is decreased from 0.025 mm 
to 0.0125 mm, however, the residual stresses do converge (Figs. 3.9 and 3.11). 
So it was decided to use the increment size of 0.025 mm in the rest of the 
study. 

It is observed that, in all the cases, the equivalent stress distribution 
(both after 20 % reduction and after unloading) in x-z plane (not shown here) 
is exactly identical to that in y-z plane. This is because the width of the block 
(dimension in y direction) is equal to its thickness (dimension in x direction). 
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y 
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Figure 3.10: Residual equivalent stress distribution after unloading m(a) x_y 
(b) yjz planes for 10 x 10 x 12 elements and increment size of 0.025 mm 
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Figure 3.11: Residual equivalent stress distribution after unloading in (a) x_y 
(b) y_z planes for 5 x 5 x 12 elements and increment size of 0.0125 mm ( \n V\pQ.’) 


3.2 Validation 


The FEM code is validated by comparing the results on forging load with 
experimental r^ults of [28]. The block material is A1-11(K) alunainium and has 
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the dimensions of 19.05 7??m x 19.05 mm x 9.525 mm. Lubricated condition 
at the interface is maintained in the experimental studies of J28]. The FEM 
analysis is done using only 1/8 of the block, 5 x 5 x 12 elements and increment 
size of 0.025 mm. Complete sticking at the interface is assumed in the FEM 
analysis. 



Figure 3.12: Comparison of theoretical and experimental [28] load- 
displacement curves. 

Fig. (3.12) shows the comparison of FEM and experimental load- 
displacement curves. The discrepancy between the theoretical and experimen- 
tal predictions can be explained as follows. In experimental studies, as the 
die-workpiece interface is lubricated, slip is observed over a part of the inter- 
face, whereas in this FEM modelling, full sticking is assumed over the whole 
interface. As a result, the elements on the free surface near the edge undergo 
a very large deformation. But in the experimental studies, the deformation 
is not so large due to the slip at the interface. Even though the FEM model 
neglects the interfacial slip, the energy needed to create the above mentioned 
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large deformation of the free surface elements seems to be more than what 
is dissipated due the interfacial slip. This is the reason why the theoretically 
predicted \-alues of the forging load are higher than the experimental values. 


3.3 Analysis of Forging Parameters for Typi- 
cal Case 

In this section, various forging parameters like the equivalent stress distribu- 
tion, residual stress distribution, contact pressure distribution, bulge profile 
etc are discussed for a typical block of 5 mm x 5 mm x 5 mm size and made up 
of Al-llOO aluminium. The mesh of 5 x 5 x 12 elements and the increment size 
of 0.025 mm used. The analysis is carried out up to 20% reduction in height. 

Contact pressure: 

Figure (3.13) shows the contact pressure distribution at the die- workpiece 
interface. The contact pressure is minimum at the center of the interface and 
increases in both the width-wdse and thickness-wise directions. It attains a 
maximum ralue at the free corner. The high \'alue at the free corner is the 
result of complete sticking at the interface. 

Equivalent stress after loading: 

Figure 3.5 shows the contours of equivalent stress in x-y and y-z planes. 
It show’s that the entire domain has yielded as every where aeq > 62.74 M pa, 
the initial }ield stress. In x-y plane, the equi\alent stress is maximum at the 
origin of X and y axes (i.e. at the center of block) and decreases in the 
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diagonal direction towards the free corner. In y-z plane also, the maximum 
value of equivalent stress occurs at the center. But, here stress decreases in 
the upward direction. The minimum value, however, does not occur along the 
entire interface but only at the center of the interface. 


It is observed that the contours of equivalent stress in x-z plane (not 
showm here) are identical to that in y-z plane. This is because the width of the 
block (dimension along y direction) is equal to its thickness (dimension along 
X direction) 

Bulge profile after loading: 

The bulge profile after loading is showm in Fig (3.14). The bulge is more 
along X and y axes than at the free corner. 
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Figure 3.14: bulge profile after loading 
Residual stresses : 

Figure 3.9 shows the contours of residual equivalent stress (after unload- 
ing) in x-y and y-z planes. The pattern of these distributions are different 
than those of equivalent stress after loading. Further, the values are less. In 
x-y plane, all the values are less than the initial yield stress (62.74 Mpa). The 
maximum value now' is not at the origin of x and y axes but somewhere in 
the middle of the plane. In y-z plane, along the interface, there are pockets of 
> 62.74 Mpa. The maximum value of residual equivalent stress occurs in 
this region. 

Due to geometric symmetry (a = b), the r^idual equivalent stress dis- 
tribution in x-z plane is exactly identical to that in y-z plane. That is why it 
is not show'n here. 
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Bulge profile after unloading: 


The bulge profile after unloading is shown in Fig. (3.15). During unload- 
ing, only the elastic deformation is recovered, therefore there is only a small 
change in shape and size from the bulge profile of Fig. 3.14 . 


z 



Figure 3.15: bulge profile after unloading 
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3.4 Parametric Study 


In this section, a study is made of the effects of process variables viz. the 
percent reduction in height, the geometric ratios h/a and b/a and the material 
properties on the equi%'alent stress distribution after loading and the residual 
stress distribution. 

The mesh of 5 x 5 x 12 elements and the incremental displacement of 

0. 025 mm is used in this study. 

3.4.1 Effects of Geometric Ratios (h/a and b/a) 

The effects of geometric ratios h/a (height to width ratio) and b/a (thickness 
to width ratio) are studied by carrying out the analysis for the following 3 
cases. Case 1 : h/a = 1, b/a = 1, Case 2 : h/a = 0.6, b/a = 1 and Case 3 : 
h/a = 1, b/a = 0.6 . The value of'^used is 5 mm. The analysis is carried out 
up to 20 % reduction cind the material used is Al-llOO aluminium. For the case 

1, the equivalent stress distribution after loading and the residual equivalent 
stress distribution both in x-y and y-z planes are already shown in Figures 3.5 
and 3.9 . The distributions in x-z plane are identical to those in y-z plane and 
hence are not shown. 

When the height of the block is reduced (case 2), then the stress dis- 
tribution in x-y plane doesn’t change. This is as expected. However, "when 
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the thickness b (dimension along y axes) of the block is reduced compared 
to its width a (dimension along x axes), i.e. in case 3, it was expected that 
the equivalent stress distribution in x-y plane would change. But it doesn't 
change. Since the stress distribution in x-y plane for all the three cases is 
identical, it is shown here only for the case 1 (Fig. 3.5) 

The equivalent stress distribution in y-z plane certainly changes for the 
case 2, i.e when the height of the block is reduced (see Figures 3.5 and 3.16). 
The maximum and minimum values change by a small amount. Because of the 
geometric symmetry (b = a), the distribution in x-z plane is exactly identical 
and that is why it is not shown. It was expected that when b is reduced 
compared to a (case 3), the stress distribution in y-z plane would change. But, 
comparison of Figure 3.5 and 3.17 shows that, it doesn’t happen. It was also 
expected that, now since the geometric symmetry does not e.xist (b T^a). the 
distribution in x-z plane would be different than that in y-z plane. But, this 
also doesn’t happen (see Fig. 3.17). 



Figure 3.16: Equivalent stress distribution after 20 % reduction in y-z plane 
for h/a=0.6 ( In 
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Figure 3.17: Equivalent stress distribution after 20 % reduction in (a) yjz (b) 
xjz planes for b/a=0.6 Cio MpcsJ) 


The pattern of residual equivalent stress distribution in x-y plane doesn’t 
change much for the case 2, i.e. when the height of the block is reduced (see 
Figures 3.9 and 3.18). However, the value of maximum stress becomes almost 
double the value for case 1. Further its location shifts to the end of y axes. 
For the case 3, i.e. when b is decreased compared to a, the residual stress 
distribution becomes unsymmetric with respect to x and y directions (see Fig. 
3.20). Also, the value of maximum residual stress increases. Further, it’s 
location shifts to the origin of x and y axes. Thus, it can be concluded that 
the values of residual stresses increase when any one dimension is decreased, 
the increase being more when the height is decreased. 

The residual equivalent stress distribution in y-z plane changes when the 
height of the block is reduced, i.e. for the case 2 (see Figures 3.9 and 3.19). 
Further, the value of maximum stress increases substantially even though its 
location doesn’t change. The distribution in x-z plane (not shown) is identical 
to that in y-z plane. For the case 3 (i.e. when b is reduced compared to 


65 





a), it was expected that the residual stress distribution in y-z plane would 
change. But, comparison of Figures 3.9 and 3.21 shows that it doesn’t happen. 
However, the distribution in x-z plane is different than that in y-z plane. The 
value of maximum residual stress is also slightly higher. 



X 


Figure 3. 1 8: Residual equivalent stress distribution after unloading in x-y plane 
for h/a=0.6 (.("n MPcl') 



Figure 3.19: Residual equivalent stress distribution after unloading in y-z plane 
for h/a=0.6 (_ 
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Figure 3.21: Residual equivalent stress distribution after unloading in (a) y-2 
(b) xjz planes for b/a =0.6 (^\r) 
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3.4.2 Effect of Percentage Reduction in Height 

The effect of percentage reduction (% r) is studied by analysing the problem 
for 2 different values of r : 20 % and 10 %. The material used is Al-llOO 
aluminium and the geometric ratios are h/a = 1 and b/a = 1 with a = 5 
mm. For 20 % reduction, the equivalent stress distribution after loading and 
residual equivalent stress distribution are already presented in Figure 3.5 and 
3.9. 


When the percentage reduction is reduced, the pattern of equivalent 
stress distribution doesn’t change. However, as expected, the value of maxi- 
mum equivalent stress decreases with a reduction in % r. At 20 % reduction, 
the value of maximum equivalent stress is 120 Mpa (Fig. 3.5) which reduces 
to 102 Mpa at 10 % reduction. 

As far as the pattern of residual equivalent stress distribution is con- 
cerned, it also doesn’t change much when the reduction is reduced to 10 7c 
(see Figures 3.22 - 3.23). How'ever, as expected, the value of maximum resid- 
ual equi\-alent stress decreases with a reduction in % r. The maximum value 
decreases from 65 Mpa (Fig. 3.9) at 20 % reduction to 55 Mpa at 10 % reduc- 
tion. Thus, at 10 % reduction, all the residual equivalent stress values are less 
than the yield stress (62.74 Mpa). 
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Figure 3.22: Residual equivalent stress distribution after unloading in x-y plane 
for 10 %reduction (. V^po.) 



Figure 3.23: Residual equivalent stress distribution after unloading in pj^ne 
for 10 %reduction ( in I^ pO 




3.4.3 Effect of Material Properties 


The effect of material properties is analysed by considering 2 different materials 
: Al-llOO aluminium and AISI-1019 steel. The analysis is done up to 20 % 
reduction and for the geometric ratios of h/a = 1 and b/a =1 (with a = 5 mm). 
The equivalent stress distribution after unloading and the residual equivalent 
stress distribution for the case of aluminium are already shown in Figures 3.5 
and 3.9 . For the steel, these distributions are shown in Figures 3.24 and 3.25. 
Comparison of these Figures shows that the pattern of both these distributions 
for the steel is the same as that for the aluminium. However, the maximum 
values are more for the steel. The maximum equivalent stress for aluminium is 
120 Mpa (Fig. 3.5) whereas for the steel, it is 1030 Mpa (Fig. 3.24). Similarly, 
the maximum residual equivalent stress for aluminum is 65 Mpa (Fig. 3.9), 
but for steel, it is 660 Mpa (Fig. 3.25). Just as in the case of aluminium, only 
in a small part of the y-z plane, the residual equi\-alent stress exceeds the yield 
stress (which is 478.2 Mpa for the steel). 
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Chapter 4 


Conclusions and Scope for 
Future work 


The objective of this work is to develop a 3-D large deformation elasto-plastic 
finite element code for determination of residual stresses in cold forging pro- 
cess. The updated Lagrangian formulation which is convenient for handling 
material and geometric nonlinearities is used. New incremental objective stress 
measure and logarithmic strain measure are used which can allow the use of 
large incremental displacement. Forging process is a displacement controlled 
problem. Therefore, to accelerate the convergence of the iterative scheme, arc 
length method is used in conjunction with modified Newton Raphson itera- 
tive technique. The detailed parametric study of residual stresses has been 
presented to show the effects of three process variables namely the percentage 
reduction in height, the geometric ratios and the material properties. For a 
typical case, the normal stress stress distribution at the die-workpiece interface 
and the deformation patterns after loading and unloading have been presented. 
The following conclusions can be drawn : 
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• After loading, the equivalent stress is maximum at the center of the block, 
but after unloading the maximum residual equiralent stress occurs at the 
die-workpiece interface. 

• Prom parametric study, it can be concluded that the values of residual 
stresses increase when any one dimension is decreased, the increase being 
more when the height is decreased. 

• The material properties have no effect on the pattern of either the equiv- 
alent stress distribution after loading or the residual stress distribution. 

• The 7c reduction also has no effect on the pattern of the equivalent stress 
distribution after loading as well as the residual stress distribution. The 
values of residual stress increases with 7c reduction. 

The following work may be taken up as extensions of the present study ; 

• Proper modelling of interfacial friction and slip. 

• Dependence of material behaviour on strain rate (viscoplasticity) and 
temperature. 

• Analysis and prediction of defects in forging proc^ during loading and 
due to residual stresses. 

• A preprocessor and postprocessor are a must for generation of input data 
and interpretation of results in a complicated 3-D analysis. 
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